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Abstract 


This  report  describes  the  technical  progress  under  a  two-year  research  effort  to  inves- 
tigate^the  properties  of  the  Lockheed  3-D  elliptic  grid  generation  procedure  and  improves n'f  ^ 
its  capabilities  for  treating  complex  geometric  features  characteristic  of  aircraft  configura¬ 
tions.  T’hetnitial  investigation  ldentified>two  areas  of  weakness  in  the  original  technique.^7 
BotFvrelate  to  the  fundamentally  important  capability  for  generating  curved-surface  grids 
as  a  necessary  prerequisite  for  generating  3-D  space  grids.  The^reSa- of -weakness  were:  •(*)  -  i ) 
difficulties  asspciated  with  the  limitations  of  Cartesian  coordinates  for  describing  convex 
surfaces;  and  (ft/  problems  in  controlling  the  orthogonality  of  the  grid  near  boundaries. 

The  bulk  of  the  first  year’s  research  3uas  focussed  on  the  first  weakness,  an  ^resulted  in 
development  of  a  generalized  system  of  elliptic  grid  generation  equations  that  allow  both 
the  surface  and  the  grid  to  be  described  in  terms  of  an  arbitrary  system  of  curvilinear 
coordinates.  "tk*  '  '  d  t  / '  <">  5  -■> 

The  second  year  of  the  contract  has  beeh^concerned  primarily  with  the  problem  of 
grid  orthogonality  near  boundaries.  Three  approaches  have  beetfr investigated  to  improve 
the  control  that  the  user  has  over  the  angle  at  which  grid  lines  intersect  the  boundaries^?  W  *  •  r 
Thethree  approaches  ar^(l)  refinement  of  the  original  techniques;  (2)  development  of 
techniques  for  marching  from  a  given  boundary  using  the  angle  as  an  initial  condition,* 
and  (3)  derivation  of  a  fourth-order  system  that  permits  additional  boundary  conditions 
on  the  grid  intersection  angle,  |  ^  j  ^  v;  a  •  ,  «•  r 
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Section  1 

INTRODUCTION 

This  report  describes  the  technical  progress  under  a  two-year  research  effort  to  inves¬ 
tigate  the  properties  of  a  3-D  elliptic  grid  generation  procedure  and  improve  its  capabilities 
for  treating  complex  geometric  features  characteristic  of  aircraft  configurations.  The  grid 
generation  technique  that  was  the  subject  of  the  research  was  an  outgrowth  of  earlier  Lock¬ 
heed  work  on  two-dimensional  grids  for  computing  the  flow  in  and  about  three-dimensional 
aircraft  jet  exhaust  nozzles  (Thomas  (1970, 1980),  Thomas  and  Middlecoff  (1980)].  In  that 
work,  the  aim  was  to  build  up  a  three-dimensional  grid  as  a  collection  of  two-dimensional 
grids  in  successive  cross-sectional  planes  along  the  nozzle.  A  method  was  developed  that 
was  patterned  after  the  earlier  work  of  Thompson  [Thompson, 1974]  in  which  a  2-D  grid  is 
generated  numerically  as  the  solution  to  an  elliptic  boundary  value  problem.  We  improved 
upon  this  approach  by  incorporating  a  set  of  general  grid  control  functions  into  the  dif¬ 
ferential  operator  of  the  elliptic  system  of  partial  differential  equations.  In  this  method, 
we  evaluated  the  grid  control  functions  at  the  boundaries  and  interpolated  them  into  the 
interior  of  the  region  to  provide  direct  control  over  the  behavior  of  the  interior  grid  by 
making  it  mimic  the  boundary  shape  and  grid  point  distribution  [Thomas,  1979;  Thomas 
and  Middlecoff,  1980].  Later,  we  extended  this  two-dimensional  elliptic  grid  generation 
technique  to  three  dimensions  (Thomas,  1981].  The  latter  technique  is  referred  to  herein 
as  the  Lockheed  3-D  elliptic  grid  generation  procedure,  and  has  been  the  subject  of  the 
two-year  research  effort  described  in  this  report. 

During  the  first  year  of  the  effort,  research  was  conducted  to  investigate  the  capability 
of  the  Lockheed  3-D  elliptic  grid  generation  procedure  for  handling  complex  geometrical 
features  characteristic  of  aircraft  configurations.  The  initial  investigation  identified  two 
areas  of  weakness  in  the  original  technique.  Both  relate  to  the  fundamentally  important 
capability  for  generating  curved-surface  grids  as  a  necessary  prerequisite  for  generating  3-D 
space  grids.  The  areas  of  weakness  were:  (i)  difficulties  associated  with  the  limitations  of 
Cartesian  coordinates  for  describing  convex  surfaces,  and  (ii)  problems  in  controlling  the 
orthogonality  of  the  grid  near  boundaries.  The  latter  prolems  are  more  difficult  to  deal  with 
in  the  case  of  two-dimensional  grids  on  curved  surfaces  than  for  planar  two-dimensional 
grids. 

The  bulk  of  the  first  year's  research  was  focussed  on  the  first  weakness,  and  resulted  in 
development  of  a  generalized  system  of  elliptic  grid  generation  equations  that  allow  both 
the  surface  and  the  grid  to  be  described  in  terms  of  an  arbitrary  system  of  curvilinear  coor¬ 
dinates  [Whitney  and  Thomas,  1983].  This  overcame  the  inherent  limitations  of  Cartesian 
coordinate  representations.  Some  progress  also  was  made  in  improving  the  orthogonality 
near  the  boundaries  by  the  ad  hoc  method  of  specifying  the  curvature  distribution  that 
enters  into  the  evaluation  of  the  grid  control  functions. 

As  a  result  of  the  first  year  of  effort,  it  was  concluded  that  the  method  had  good 
potential  for  applications.  In  fact,  other  related  activities  included  application  of  the 
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method  to  practical  cases,  such  as  a  3D  grid  over  a  wing-fuselage  combination,  for  use  with 
the  FLO-57  transonic  flow  code  [Thomas  and  Neier,  1983],  and  development  of  software 
for  interactive  application  of  the  technique. 

The  second  year  of  the  contract  has  been  concerned  primarily  with  the  problem  of 
grid  orthogonality  near  boundaries.  Three  approaches  have  been  investigated  to  improve 
the  control  that  the  user  has  over  the  angle  at  which  grid  lines  intersect  the  boundaries. 
The  three  approaches  are  (1)  refinement  of  the  original  techniques,  (2)  development  of 
techniques  for  marching  from  a  given  boundary  using  the  angle  as  an  initial  condition, 
and  (3)  derivation  of  a  fourth-order  system  that  permits  additional  boundary  conditions 
on  the  grid  intersection  angle. 

A  study  has  been  made  to  examine  various  ways  of  refining  the  original  grid  control 
technique  to  allow  direct  control  over  the  angle  with  which  the  transverse  grid  lines  meet 
the  boundaries  of  the  region  in  which  the  grid  is  generated.  These  are  the  boundaries 
at  which  the  grid  control  parameters  in  the  elliptic  system  are  evaluated  in  terms  of  a 
pre-assigned  distribution  of  boundary  points.  One  effective  refinement  has  been  developed 
that  is  based  on  the  fact  that  the  equation  defining  each  grid  control  parameter  at  a 
boundary  is  composed  of  two  terms,  a  grid  curvature  term  and  a  grid  stretching  term. 
In  the  refined  technique,  the  two  terms  are  evaluated  on  different  parts  of  the  boundary, 
are  interpolated  separately  into  the  interior  of  the  region,  and  are  summed  to  form  the 
grid  control  parameter.  In  contrast,  the  original  technique  evaluated  both  terms  along  the 
same  part  of  the  boundary,  and  the  grid  control  parameter  itself  was  interpolated  into  the 
interior.  The  refined  technique  improves  considerably  the  user's  control  of  the  grid,  but 
also  has  some  limitations:  (1)  it  does  not  change  the  local  behavior  of  interior  grid  lines 
when  the  boundary  segments  are  straight,  (2)  it  ensures  orthogonality  of  grid  lines  to  a 
boundary  only  when  the  two  boundary  segments  belonging  to  the  same  family  as  the  grid 
lines  are  themselves  both  curved  and  orthogonal  to  the  boundary  in  question,  and  (3)  the 
technique  is  less  useful  for  two-dimensional  grids  on  curved  surfaces  than  for  either  planar 
two-dimensional  grids  or  three-dimensional  space  grids. 

The  marching  approach  is  based  on  using  our  original  second-order  system  of  elliptic 
equations,  including  the  grid  control  functions,  and  solving  the  equations  by  forward  in¬ 
tegration  away  from  an  initial  boundary  at  which  both  the  grid  point  locations  and  grid 
intersection  angles  are  specified  as  initial  conditions.  This  approach  was  abandoned  when 
we  were  unable  to  develop  a  stable  marching  procedure  for  general  grid  topologies. 

As  discussed  previously,  the  original  second-order  system  cannot  accommodate  speci¬ 
fied  grid  intersection  angles  at  the  boundaries.  To  overcome  this  shortcoming,  one  approach 
is  to  go  to  a  higher-order  system  of  partial  differential  equations  which  admit  more  bound¬ 
ary  conditions.  We  have  investigated  such  a  fourth-order  system  obtained  by  applying 
the  Laplacian  differential  operator  to  the  original  second-order  system.  This  particular 
system  has  the  advantage  of  retaining  the  grid  control  functions  inherent  in  the  original 
second-order  system  while  allowing  the  grid  intersection  angles  at  boundaries  to  be  speci¬ 
fied  as  boundary  conditions.  The  resulting  system  is  significantly  more  complex  than  the 
second-order  system.  Numerical  results  have  been  obtained  for  several  representative  cases 
in  two  dimensions,  but  further  effort  would  be  required  to  develop  a  reliable  and  efficient 
algorithm  for  solving  the  fourth-order  equations  in  three  dimensions. 
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The  remainder  of  this  report  contains  a  detailed  discussion  of  the  described  research. 
The  objectives  of  each  year  of  the  effort  are  given  in  Section  2  as  they  appeared  in  the 
original  contract  statements  of  work.  The  research  performed  toward  those  objectives  is 
discussed  in  Section  3.  Finally,  Section  4  lists  the  technical  papers  prepared  under  the 
contract  and  the  personnel  who  participated  in  the  research. 


Section  2 

OBJECTIVES  OF  THE  RESEARCH 


The  research  objectives  of  each  year  of  the  research  effort  are  given  below  as  they 
appeared  in  the  contractual  statements  of  work. 

2.1  Objectives  of  the  First  Year  of  Effort 

Investigate  the  capability  of  the  Lockheed  3-D  elliptic  grid  generation  procedure  in 
handling  complex  geometrical  features  which  are  characteristic  of  aircraft  configurations. 
Determine  limitations  in  the  procedure  by  considering  model  problems  which  include  spe¬ 
cific  geometric  features  which  may  violate  known  or  suspected  restrictions  in  the  method. 
Eliminate  current  constraints  which  exclude  non-orthogonally  intersecting  surfaces,  such 
as  the  intersection  of  a  swept  wing  with  a  fuselage. 

2.2  Objectives  of  the  Second  Year  of  Effort 

Investigate  ways  to  control  the  angle  at  which  grid  lines  intersect  boundaries.  Three 
approaches  will  be  considered:  (1)  refinement  of  grid  control  techniques  on  the  current 
method,  (2)  development  of  march  solution  techniques  that  would  allow  the  angle  of  inter¬ 
section  with  the  boundary  to  be  prescribed  as  an  initial  condition,  and  (3)  derivation  of  a 
fourth-order  elliptic  system  based  on  the  present  second-order  system  to  provide  additional 
boundary  conditions  on  the  grid  angle  at  the  surface. 
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Section  3 

SUMMARY  OF  THE  RESEARCH  ACCOMPLISHED 

The  following  paragraphs  summarize  the  research  performed  in  each  year  of  the  effort 
toward  the  objectives  stated  in  Section  2. 

3.1  Research  Performed  During  the  First  Year  of  Effort 

The  effort  during  the  first  year  was  focussed  primarily  on  the  use  of  elliptic  equation 
systems  to  generate  quasi-two-dimensional  grids  on  curved  surfaces  typical  of  the  surfaces 
found  on  aircraft  configurations.  The  ability  to  generate  smooth,  well-controlled  surface 
grids  is  a  necessary  prerequisite  for  three-dimensional  grid  generation  by  elliptic  differential 
systems,  because  the  surface  grid  on  the  boundary  of  a  3-D  region^ forms  the  boundary 
values  that  are  needed  to  solve  the  3-D  elliptic  equation  system  that  generates  the  3-D 
grid. 

Earlier,  we  had  developed  a  special  elliptic  equation  system  for  generating  grids  on 
surfaces  z=f(x,y)  in  Cartesian  coordinates  [Thomas,198l].  This  system  is  similar  to  the  one 
we  had  used  previously  for  plane  2-D  grids  [Thomas  and  Middlecoff  (1979)],  but  contains 
forcing  function  "source  terms"  that  involve  the  slope  and  curvature  of  the  surface.  During 
the  first  year,  two  difficulties  were  discovered  in  the  use  of  this  Cartesian-based  elliptic 
system:  (i)  problems  in  obtaining  a  converged  solution  for  the  grid  on  fuselage-like  or  wing¬ 
like  surfaces  in  regions  where  the  derivatives  of  the  surface  shape  function  f(x,y)  are  large  or 
are  rapidly-varying,  and  (ii)  difficulties  with  the  local  control  of  grid  lines  near  boundaries, 
where  near-orthogonality  of  the  grid  lines  to  the  boundary  usually  is  desirable,  and  often  is 
absolutely  necessary  to  assure  smoothness  of  the  grid  across  boundaries  between  different 
regions. 

The  first  problem  is  a  consequence  of  using  Cartesian  coordinates  in  terms  of  which  the 
surface  shape  function  f(x,y)  is  always  multiple-valued  and  has  regions  of  locally  unbounded 
derivatives  when  used  to  describe  a  convex  geometrical  object  such  as  an  aircraft  or  its 
components  (fuselage,  wings,  nacelles,  etc).  It  became  apparent  to  us  that  this  problem 
could  be  avoided  by  using  other  coordinate  systems,  such  as  cylindrical,  spherical,  etc., 
that  are  both  more  natural  and  more  convenient. 

The  major  part  of  the  first  year’s  effort,  therefore,  was  applied  to  develop  a  generalized 
system  of  elliptic  surface  grid  generation  equations  that  allows  one  to  describe  the  surface 
and  generate  the  grid  in  terms  of  any  arbitrary  curvilinear  coordinate  system.  Our  earlier 
Cartesian-based  surface  grid  generation  system  is  simply  a  special  case  of  the  generalized 
system.  The  results  of  this  work  were  reported  in  a  technical  paper  at  the  Mini-Symposium 
on  Advances  in  Grid  Generation  [Whitney  and  Thomas,  1983]. 

The  second  problem  mentioned  above,  namely  grid  control,  also  was  addressed  initially 
during  the  first  year.  It  was  found  that  our  original  techniques  of  controlling  the  behavior 
of  the  grid  through  the  boundary  values  for  the  elliptic  system  are  much  less  powerful 
for  dealing  with  grids  on  curved  surfaces  than  we  had  found  them  to  be  for  plane  grids. 
An  example  is  given  in  Fig.  1  to  demonstrate  the  inability  of  the  original -*gr  id  control 
technique  to  produce  grid  lines  that  are  orthogonal  to  the  boundaries.  The  figure  shows  a 
grid  on  the  surface  of  a  bilaterally  symmetric  fuselage,  a  spherically-capped  cylinder  with  a 
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cutout  for  a  mid-wing  having  a  NACA  0012  profile.  The  grid  was  generated  by  the  elliptic 
equation  system  mentioned  earlier,  using  the  grid  control  techniques  originally  developed 
for  plane  grids.  One  can  see  that  the  grid  lines  are  quite  skewed  and  do  not  intersect  the 
longitudinal  boundaries  orthogonally  except  very  near  the  nose  tip  and  again  very  near  the 
aft  end.  Clearly,  if  this  grid  for  the  upper  half  of  the  fuselage  were  joined  to  a  similar  one 
for  the  lower  half,  the  "vertical”  grid  lines  would  have  a  sharp  discontinuity  in  slope  as  well 
as  in  curvature  in  crossing  the  boundary  where  the  two  halves  join.  A  similar  skewness 
occurs  along  the  boundary  between  the  right  and  left  halves,  the  former  of  which  is  not 
shown. 

Some  progress  was  made  during  the  first  year  to  improve  the  orthogonality  near  the 
boundaries  by  specifying  the  curvature  distribution  that  enters  into  one  of  the  two  terms 
that  are  used  in  evaluating  the  grid  control  functions  at  each  boundary.  The  grid  shown 
in  Fig.  2  was  produced  in  this  way,  and  has  transverse  grid  lines  that  are  very  nearly 
orthogonal  to  all  boundaries,  including  those  where  the  original  grid  of  Fig.  1  is  markedly 
skewed. 

3.2  Research  Performed  During  the  Second  Year  of  Effort 

Three  different  approaches  were  investigated  during  the  second  year  to  improve  control 
over  the  angle  with  which  the  grid  lines  intersect  boundaries:  refinement  of  the  existing  grid 
control  technique,  marching  methods  that  would  allow  the  angle  of  intersection  with  an 
initial  boundary  to  be  prescribed  as  an  initial  condition,  and  manipulation  of  the  original 
second-order  system  to  derive  a  fourth-order  elliptic  system  that  would  permit  additional 
boundary  conditions  on  the  intersection  angle  at  all  boundaries.  Progress  in  each  of  these 
areas  is  outlined  below. 

3.2.1  Refinement  of  the  Original  Technique 

We  have  studied  various  ways  of  refining  the  original  grid  control  technique  to  improve 
control  over  the  angle  with  which  the  transverse  grid  lines  meet  the  boundaries  of  the 
region  in  which  the  grid  is  generated.  These  are  the  boundaries  at  which  the  grid  control 
parameters  in  the  elliptic  system  are  evaluated  in  terms  of  a  pre-assigned  boundary  grid 
point  distribution.  To  simplify  the  discussion,  we  shall  use  a  two-dimensional  example  to 
contrast  the  original  technique  with  the  newly-developed  refinements. 

Overview  of  the  Original  Technique.  Given  a  simply-connected  region  of  the  x,y 
plane  such  as  the  region  illustrated  in  Fig.  3,  the  original  technique  generates  a  boundary- 
conforming  curvilinear  coordinate  transformation 

x,y  * — ♦ 

£  =  £(*,v),  *?  =  n(z,v) 

that  maps  the  region  onto  a  rectangle  (Fig.  4).  The  transformation  is  generated 
numerically  as  the  solution  to  the  elliptic  system  of  equations 


vse  =  *(e,*)  ivei’. 


v*n  =  |v»i! 
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where  <t>,  rp  are  grid  control  functions.  The  roles  of  dependent  and  independent  vari¬ 
ables  in  these  equations  are  interchanged  to  obtain  an  elliptic  system  that  can  be  solved 
numerically  on  a  uniform  rectangular  grid  in  the  domain  of  Fig.  4  to  obtain  the  grid  point 
coordinates  in  the  physical  region 

a[r(t  +  -  20F(ti  +  Tf(rnn  +  0r„)  =  0  (2) 

r  =  (x,y)  (3a) 

« =  l^l2.  '•»)»  'I  =  l^-fl2  (36) 

Boundary  values  fy*  are  assigned  along  the  perimeter  O’A’B’C’O’  of  the  square.  These 
boundary  values  represent  the  coordinates  r  =  (z,  y)  of  grid  points  along  the  perimeter 
OABCO  of  the  physical  region. 

In  the  original  grid  control  technique,  limiting  forms  of  Eq.  (2)  were  used  to  evaluate 
the  control  functions  <f>,ip  locally  along  the  boundary  of  the  square,  so  that  the  elliptic 
system  then  could  be  solved  numerically.  The  functions  were  evaluated  individually  along 
«ach  of  the  four  segments  of  the  boundary  as  follows. 

Consider,  for  example,  the  lower  horizontal  segment  O’A’,  which  represents  the  bound¬ 
ary  OA  in  Fig.  3,  and  is  a  coordinate  curve  17  =  constant  along  which  (  alone  varies.  A 
unique  equation  for  evaluating  <f>  on  this  segment  follows  by  taking  the  projection  of  Eq. 
(2)  onto  the  vector  r(  that  is  tangent  to  the  segment,  and  invoking  the  orthogonality 
constraint 

r€r,  =  0  on  O'A’B'C’O '  (4) 

The  resulting  equation  can  be  cast  in  the  form 


4  —  ~  (To  +  Ir^jTi)  (5a) 

To  =  ft  •  F{e)/|r-{|2  =  (|f{f  ){/(2|f{|2)  =  (!n|fel)f  (56) 

r,  =  (r4-r„)/(|f£||f,|J)  (5e) 

The  term  7o  is  simply  the  logarithmic  derivative  of  arc  length  with  respect  to  the 
parameter  (  along  the  boundary  segment  [Thomas,  1981].  It  depends  on  the  spacing 
between  grid  points  along  the  boundary  segment  OA,  and  will  be  referred  to  as  the  "grid 
stretching”  term.  Both  of  the  (-differentiated  quantities  |f(|  and  To  in  Eqs:’  (5)  can  be 
evaluated  numerically  from  the  boundary  values  assigned  on  the  segment  O’ A’.  However, 
the  term  T\  involves  rj -derivatives  which  cannot  be  evaluated  locally  using  only  boundary 


data  on  O’A*.  This  term  represents  the  local  curvature  (Thomas,  1981]  of  the  grid  lines 
£  =  constant  that  are  transverse  to  the  boundary  segment  OA  and  will  be  referred  to  as 
the  "grid  curvature"  term. 

Clearly,  one  could  specify  any  desired  curvature  distribution  function  Tj(£)  along 
the  boundary  segment.  We  shall  return  to  this  idea  later.  In  the  original  grid  control 
technique,  however,  this  curvature  term  was  simply  evaluated  at  the  end  points  of  the 
segment  OA,  and  was  given  a  linear  variation  over  the  rest  of  the  segment  for  use  in 
evaluating  the  grid  control  function  4>  from  Eq.  (5a).  The  same  procedure  was  employed 
to  evaluate  the  parameter  ^(£)  along  the  upper  boundary  segment  B’C’  in  Fig.  4.  A 
continuous  representation  ^(£,q)  throughout  the  interior  of  the  computational  square  was 
then  obtained  by  linear  interpolation  as  a  function  of  f)  along  vertical  mesh  lines  £  = 
constant. 

Outline  of  The  Refined  Technique.  Refinements  have  been  added  to  the  original 
technique  to  improve  the  grid  control  both  in  the  interior  of  the  region  and  near  the 
boundaries  when  the  boundaries  are  curved.  The  refinements  are  based  on  a  new  approach 
that  makes  direct  use  of  the  physical  interpretation  of  the  term  7\  in  Eq.  (5b)  as  the 
curvature  of  r; -directed  grid  lines.  The  rationale  is  as  follows. 

The  purpose  of  the  grid  control  functions  is  to  make  the  elliptic  system  (2) 
produce  an  interior  grid  that  mimics  the  shape  and  curvature  of  the  boundaries  and  the 
grid  point  distribution  that  is  assigned  on  those  boundaries.  Consider  a  point  P’  at  which 
the  grid  control  function  is  to  be  evaluated  (see  Fig.  4),  and  let  P  be  the  image  of  that 
point  in  the  physical  region  of  Fig.  3.  The  horizontal  and  vertical  grid  lines  in  Fig.  4  that 
pass  through  the  point  P’  have  as  their  image  in  the  physical  region  two  intersecting  grid 
lines  that  are  illustrated  by  the  dashed  lines  in  Fig.  3.  If  these  two  physical  grid  lines 
were  known  explicitly,  then  one  could  evaluate  the  function  <f>  locally  from  Eq.  (5).  The 
first  term  could  be  calculated  if  one  knew  only  the  £ -directed  grid  line  through  P;  this  is 
the  grid  stretching  term  that  reflects  the  grid  spacing  along  the  £-directed  grid  line.  The 
second  term  represents  the  curvature  of  the  17-directed  grid  line  through  P,  but  requires 
knowledge  of  both  grid  lines  since  it  contains  derivatives  with  respect  to  both  £  and  q. 
However,  one  can  make  use  of  the  orthogonality  constraint  (4)  to  rewrite  Eq.  (5c)  in  terms 
of  t}  derivatives  alone,  as 


(6) 


where  r  =  (x,y,z)  with  z  =  0  and  k  =  (0,0,1)  is  the  Cartesian  unit  vector  in  the 
z-direction  normal  to  the  x,  y  plane.  Eq.  (6)  follows  from  the  fact  that  the  orthogonality 
relation  (4)  can  be  written  in  the  form 


IL  =  ^lx* 

ftl  I'*! 


(7) 


The  foregoing  conceptual  exercise  suggests  a  rational  way  to  compute  the  grid  control 
function  ^  when  the  grid  is  known  only  on  the  boundaries,  and  one  would  like  the  interior 


9 


grid  produced  by  the  elliptic  system  to  be  some  smooth  interpolation  from  the  known  grid 
on  the  boundaries:  the  grid-stretching  term  To  governs  the  local  grid  spacing  along  the 
(-directed  grid  line  through  the  point  P,  and  should  be  estimated  in  terms  of  the  known 
grid  spacing  on  the  (-directed  boundary  segments  OA  and  CB.  The  curvature  term  Ti 
governs  the  local  curvature  of  the  ^-directed  grid  line  through  P  and  should  be  estimated 
in  terms  of  the  known  curvature  of  the  ^-directed  boundary  segments  OC  and  AB. 

The  refined  grid  control  procedure  then  is  as  follows.  We  compute  the  grid  control 
function  4>  at  a  general  interior  point  P’  of  the  square  in  Fig.  4  by  using  Eq.  (5a)  locally. 
The  grid  stretching  term  To  is  computed  at  the  boundaries  O’A’  and  C’B’  and  interpolated 
linearly  as  a  function  of  r?  into  the  interior  of  the  square.  The  factor  |r*f|  in  Eq.  (5a)  is 
evaluated  in  the  same  way.  The  curvature  term  is  computed  at  the  boundaries  OC  and 
AB  by  using  Eq.  (6),  and  is  interpolated  linearly  as  a  function  of  (  into  the  interior  of  the 
square. 

The  described  refinement  greatly  improves  both  the  behavior  of  the  interior  grid  and 
the  orthogonality  of  grid  lines  at  the  boundaries  whenever  the  boundary  segments  are 
curved  and  are  nearly  orthogonal  to  one  another.  The  improvement  is  illustrated  by  the 
example  shown  in  Fig.  5  for  a  region  that  has  two  curved  boundary  segments  that  are 
nearly  orthogonal  to  the  horizontal  straight  segment.  One  sees  from  Fig.  5a  that  the 
original  technique  produces  interior  grid  lines  that  roughly  follow  the  shapes  of  the  two 
curved  boundary  segments,  but  that  have  too  little  curvature  and  do  not  mimic  the  orthog¬ 
onality  with  which  the  curved  boundary  segments  meet  the  horizontal  straight  segment. 
In  contrast,  the  curved  interior  grid  lines  produced  by  the  refined  grid  control  technique 
follow  faithfully  the  curvature  of  the  boundary  segments,  and  also  are  nearly  orthogonal  to 
the  horizontal  straight  segment.  This  has  the  added  effect  of  making  the  spacing  between 
adjacent  interior  curved  grid  lines  follow  closely  the  spacing  of  the  assigned  grid  points 
along  the  straight  boundary  segments. 

The  refined  grid  control  technique  in  general  gives  similar  improvements  both  in  two- 
dimensional  plane  grids  and  in  three-  dimensional  space  grids  whenever  the  boundaries  are 
curved.  However,  there  are  some  limitations. 

1.  The  refinements  do  not  change  the  behavior  of  interior  grid  lines  when  the  boundary 
segments  are  straight.  For  example,  there  is  little  difference  between  Figs.  5a  and  5b  in 
the  degree  of  orthogonality  with  which  the  "radial”  interior  grid  lines  meet  the  upper 
(right-hand)  curved  boundary,  because  the  "radial”  boundary  segments  (the  horizontal 
and  vertical  boundary  segments)  have  no  curvature. 

2.  The  refined  technique  cannot  produce  a  family  of  interior  grid  lines  that  are  or¬ 
thogonal  to  a  boundary  unless  the  two  boundary  segments  in  that  family  are  both  curved 
and  orthogonal  to  the  boundary  in  question. 

3.  The  refined  technique  is  useful  (within  the  foregoing  two  limitations)  for  either 
planar  two-dimensional  grids  or  three-dimensional  space  grids,  but  offers  less  generality 
in  improving  two-dimensional  curved-surface  grids  generated  by  the  2-D  elliptic  systems 
developed  the  first  year  of  the  contract  (Whitney  and  Thomas,  1983]. 

The  loss  of  generality  in  the  case  of  2-D  curved  surface  grids  stems  from  the  fact 
that  the  curvature  term  T\  in  Eq.(5c)  does  not  embody  the  actual  principal  curvature  of 
the  space  curve  that  bounds  a  general  2-D  curved  surface  region  as  it  does  for  the  planar 


surface.  T\  actually  represents  only  the  component  of  the  curvature  vector  that  lies  in  the 
tangent  plane  to  the  surface,  and  is  generally  equal  to  the  principal  curvature  only  in  the 
planar  case.  This  means  that  the  effectiveness  of  the  grid  control  functions  ip  depends  on 
the  detailed  nature  of  the  curved  surface  and  on  the  shape  of  the  space  curve  that  bounds 
the  part  of  the  surface  on  which  a  grid  is  to  be  generated. 

As  discussed  in  Section  3.1,  we  have  been  able  to  enhance  the  grid  control  for  some 
specific  cases  of  surface  geometries  and  grid  topologies  by  specifying  a  priori  the  curvature 
distribution  such  as  7i(£)  along  the  f)  =  constant  lines  in  the  computational  domain  of 
Fig.  4.  We  have  found  that  this  introduction  of  ad  hoc  curvature  control  can  be  applied 
on  a  trial  and  error  basis  to  force  the  grid  to  be  nearly  orthogonal  at  boundaries,  at  least 
for  relatively  simple  surface  geometries  and  grid  topologies.  However,  the  ad  hoc  approach 
lacks  the  generality  of  the  marching  and  fourth-order  system  approaches  discussed  below. 

3.2.2  Marching  Approach 

The  marching  approach  is  closely  related  with  the  use  of  parabolic  equations  to  con¬ 
struct  the  grid  by  solving  an  initial  value  problem  [Nakamura  (1982)].  The  parabolic 
marching  idea  is  being  pursued  by  other  investigators,  and  is  largely  concerned  with  devis¬ 
ing  a  useful  system  of  equations  and  with  learning  how  to  control  the  behavior  of  the  grid 
everywhere,  not  just  near  boundaries.  In  contrast,  the  marching  approach  investigated 
under  the  contract  was  based  on  our  original  second-order  system  of  elliptic  equations, 
including  the  grid  control  functions  that  are  evaluated  at  all  boundaries.  Because  the  grid 
generation  equations  are  of  second  order,  one  could,  in  principle,  specify  two  boundary 
conditions  at  the  initial  boundary  from  which  the  solution  is  to  be  marched:  one  boundary 
condition  on  the  grid  point  locations  along  the  boundary,  and  a  second  on  the  angle  at 
which  the  grid  lines  intersect  the  boundary.  This  would  allow  complete  control  over  the 
grid  behavior  at  the  initial  boundary  by  sacrificing  a  degree  of  control  over  the  behavior 
at  the  far  boundary  toward  which  the  solution  is  marched. 

Marching  solutions  are  rarely  attempted  for  elliptic  equations  because  the  initial¬ 
boundary  value  problem  is  ill-posed  and  numerical  stability  is  difficult  to  maintain.  In 
our  original  elliptic  system,  however,  global  information  about  the  boundary  values  is 
contained  in  the  grid  control  functions  that  appear  in  the  equations  themselves.  It  seemed 
reasonable  to  hypothesize  that  this  global  boundary  information  content  might  alter  the 
character  of  the  equations  enough  to  allow  stable  solutions  by  marching  away  from  an 
initial  boundary.  Unfortunately,  this  hypothesis  did  not  prove  true  in  practice.  Numerical 
experiments  were  conducted  in  which  a  locally  linearized  implicit  marching  technique  was 
attempted  for  planar  2-D  grids.  We  were  unable  to  obtain  stable  marching  solutions  for 
general  grid  topologies,  and  abandoned  the  marching  approach  in  favor  of  the  fourth  order 
system  approach  discussed  next. 

3.2.3  Fourth-Order  System  Approach 

The  idea  of  using  a  fourth-order  elliptic  grid  generation  system  that  allows  the  freedom 
to  specify  boundary  conditions  on  the  grid  point  locations  and  on  the  angle  of  intersection 
of  grid  lines  with  the  boundary  has  been  examined  to  a  limited  extent  [Shubin,  et  al. 
(1982)].  That  study  considered  only  the  use  of  the  biharmonic  equation,  a  linear  equation 
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that  satisfies  no  maximum  principle.  The  investigators  noted  that  unacceptable  singular 
transformations  with  grid  lines  thut  cross  others  of  the  same  family  were  encountered  even 
in  relatively  simply-shaped  regions.  It  is  likely  that  this  difficulty  would  be  much  more 
probable  for  more  complex  geometries  typical  of  aircraft,  and  would  be  much  harder  to 
overcome. 

Our  original  second-order  nonlinear  elliptic  systems  obey  a  maximum  principle  that 
ensures  a  non-singular  grid.  We  have  investigated  ways  to  develop  a  fourth-order  system 
with  this  property  by  applying  a  further  second-order  linear  differential  operator  to  the 
existing  second-order  nonlinear  elliptic  system.  It  is  plausible  to  conjecture  that  the  exis¬ 
tence  of  a  maximum  principle  for  the  present  basic  second-order  system  should  imply  that 
a  fourth-order  system  derived  from  it  also  would  obey  the  same  maximum  principle.  Such 
a  fourth-order  system  would  retain  the  grid  control  functions  and  therefore  the  quality  of 
interior  grid  control  that  is  obtained  with  the  second-order  system,  and  would  also  allow 
complete  local  control  over  the  angle  at  which  grid  lines  intersect  the  boundary. 

The  system  that  immediately  suggests  itself  is  obtained  by  operating  on  the  original 
second-order  system  with  the  Laplacian  operator;  i.eM 

V2{V2{,-*,  |V  £,|2}  =  0,  i=  1,2,3  (8) 

where  (£i»  62,^3)  —  (£>*?,?)  are  the  transformed  coordinates,  and  = 

(^,V»,  w)  are  the  grid  control  parameters.  Just  as  for  the  second-order  system,  the 
transformed  syBtem  of  equations  is  derived  by  interchanging  the  roles  of  dependent  and 
independent  variables.  Rather  than  working  with  a  single  fourth-order  system  for  each 
grid  variable,  it  is  convenient  to  consider  the  following  coupled  second-order  system: 

V2fc  =  *,  |Vfc|2  +  *,  V2p,  =  0,  .  =  1,2,3  (9) 

The  functions  p,  are  dummy  variables  that  serve  to  reduce  the  order  of  the  original 
system  in  Eq.  (8).  It  is  easily  verified  that  Eq.  (8)  is  obthained  by  eliminating  p,  from 
Eq.  (9).  The  system  in  Eq.  (9)  transforms  to  the  quasi-linear  system 


•<> 


in  which  the  coefficients  and  the  Jacobian  J  are  given  by 


In  the  definition  of  Cu  ,  the  indices  (t,j,k)  and  (/, m, n)  occur  in  cyclic  order.  These 
equations  hold  for  each  physical-space  coordinate,  r  =  (z,y,z)  ,  and  for  each  "dummy 
variable”,  p  =  (p],P2>ls)*  There  are,  therefore,  a  total  of  six  scalar  equations.  These 
equations  bear  a  close  resemblance  to  the  original  equation  system,  which  may  be  written 

as 


_  fd2f  dr  \  _  _  d2f 

£c<<(ae?  +  +  2  E  ■ 

»<i 


(12) 


The  only  difference  is  the  occurrence  of  the  J 2  term  in  Eqs.  (10)  and  (11).  With 
slight  modification,  the  same  SLOR  solution  technique  used  on  the  oritinal  system  (12) 
may  be  employed  here,  only  now  there  are  are  six  second-order  equations  instead  of  three. 

It  remains  to  determine  appropriate  boundary  conditions.  Boundary  conditions  for 
r  are  the  same  as  for  the  original  second  order  system  (12);  i.e.,  f  is  specified  at  each 
boundary  node.  Grid  intersection  slopes  on  the  boundaries  may  be  specified  by  also  giving 
the  normal  derivative  of  f  at  the  boundary.  However,  since  the  system  Eqs.  (10)  and 
(11)  is  second  order,  only  boundary  conditions  for  r  and  p  need  be  given.  A  Dirichlet 
boundary  condition  for  each  component  of  p  can  be  found  from  the  first  equation  in  Eq. 
(9),  which  transforms  to  Eq.  (10).  The  boundary  condition  is  applied  in  the  context  of  an 
SLOR  iteration  loop,  for  which  it  is  assumed  that  an  approximate,  or  guessed,  value  of  r 
is  known  at  each  node  for  a  particular  iteration  step.  The  boundary  condition  is  obtained 
from  the  limiting  form  of  Eq.  (10)  at  the  boundaries.  For  example,  consider  the  boundary 
condition  on  £3  =  f  =  1.  On  this  surface  suppose  that  r  and  rf  are  given  input  data.  Then 
by  differentiation,  all  the  other  first  and  second  partial  derivatives  of  r  are  known  except 
for  fjf.  But  this  can  be  estimated  from  r  and  fs  at  the  boundary  plus  r  on  the  first  layer 
of  nodes  interior  to  the  boundary;  i.e., 


(13) 


With  r*ff  so  approximated,  and  similar  expressions  on  the  other  boundaries,  the  three 
scalar  equations  represented  by  Eq.  (10)  may  be  used  to  solve  for  the  three  p,’i s  at  the 
boundaries.  Interior  values  of  p  are  found  by  solving  Eq.  (11). 

A  numerical  solution  to  the  fourth-order  system  of  equations  has  been  obtained  for 
several  cases.  The  numerical  method  employs  finite-difference  appraximants  of  Eqs.  (10) 
and  (11)  that  incorporate  Dirichlet  boundary  conditions  for  r  and  p  .  The  solution  then  is 
obtained  by  an  iterative  technique  based  upon  the  successive  line  over-relaxation  (SLOR) 
algorithm.  Numerical  results  have  been  obtained  for  several  two-dimensional  cases,  but 
further  effort  would  be  required  to  develop  a  reliable  algorithm  for  solving  the  fourth- 
order  system  in  three  dimensions.  The  two-dimensional  numerical  results  are  given  in 
(Whitney,  1985],  along  with  a  complete  derivation  of  the  fourth-order  system  and  a  detailed 
description  of  the  iterative  solution  technique. 
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Figure  4  Rectangular  Computational  Domain 
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Comparison  of  the  original  and  the  new  refined  grid  control  techniques 
for  a  planar  region  with  curved  boundaries. 
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